//Ghilardi and Zilberman (2021)

// Endogenous Variables

var C p Div DivTau Div_Y Div_K W Rk K N I q Y TauK TauD TauI theta A phi q_GDP q_GDP_a p_GDP p_GDP_a YY CC II KK qq pp DivDiv DivTau_annual Div_Y_annual Div_K_annual IK IK_a IY IY_a CY CY_a AA theta_annual TauD_annual phi_annual N_annual;

// Exogenous Variables (Shocks)

varexo epsA epstheta epsTauD ;

// Parameters and Calibration

parameters betta delta alp N_ss theta_ss TauD_ss TauN_ss TauI_ss TauK_ss rhoA rhoTauD rhotheta A_ss gamma sd_theta sd_A sd_TauD 
phi_ss q_ss K_ss I_ss Y_ss C_ss h_ss W_ss Rk_ss Div_ss DivTau_ss p_ss I_y_ss K_y_ss I_k_ss C_y_ss Div_y_ss DivYtau Div_K_ss q_GDP_ss p_GDP_ss IY_ss CY_ss;
 
betta=0.96;
delta=0.083;
alp=0.30;
N_ss=0.3;
theta_ss=0.0914;
TauD_ss=0.16;
TauN_ss=0; 
TauI_ss=0.0549;
TauK_ss=0.35;
rhoA=0.90; 
rhoTauD=0.99; 
rhotheta=0.99;
A_ss=1;
gamma=0.54;
sd_theta=0.14;
sd_A=1;
sd_TauD=1;

//Steady State Model

/STEADY STATE LAGRANGE MULTIPLIER ON BORROWING CONSTRAINT

phi_ss=(delta/theta_ss)-(1-TauD_ss)*(1+TauI_ss);       //Use when considering a binding steady-state
//phi_ss=0;                                            //Use when the credit constraint is always slack. 

//REST OF STEADY STATE VALUES

q_ss=phi_ss+(1-TauD_ss)*(1+TauI_ss);
K_ss=(((alp*(1-TauK_ss))/((1+TauI_ss+((phi_ss)/(1-TauD_ss)))*((1/betta)-1)+delta*(1+TauI_ss)))^(1/(1-alp)))*N_ss;
I_ss=delta*K_ss;
Y_ss=A_ss*(K_ss)^(alp)*(N_ss)^(1-alp);
C_ss=Y_ss-I_ss;
h_ss=(1/C_ss)*(1-TauN_ss)*(1-alp)*(Y_ss/N_ss);
W_ss=(1-alp)*A_ss*((K_ss/N_ss)^alp);
Rk_ss=alp*A_ss*(N_ss/K_ss)^(1-alp);
Div_ss=(1-TauK_ss)*(Y_ss-W_ss*N_ss)-(1+TauI_ss)*I_ss;
DivTau_ss=(1-TauD_ss)*Div_ss;
p_ss=(betta*(1-TauD_ss)*Div_ss)/(1-betta);
I_y_ss=I_ss/Y_ss;
K_y_ss=K_ss/Y_ss;
I_k_ss=I_ss/K_ss;
C_y_ss=C_ss/Y_ss;
Div_y_ss=Div_ss/Y_ss;
DivYtau=(1-TauD_ss)*Div_y_ss;
Div_K_ss=DivTau_ss/K_ss;
q_GDP_ss=q_ss/Y_ss;
p_GDP_ss=p_ss/Y_ss;
IY_ss=I_ss/Y_ss;
CY_ss=C_ss/Y_ss;

model;

//NON-BINDING MODEL - Use these 2 equations when simulating the permanently slack model with phi_ss=0

//phi=0;

//q=(1-TauD)*(1+TauI+gamma*((I/K(-1))-delta));


//OBC MODEL - Use these 2 equations when allowing for OBC and phi_ss>0.  

phi=q-(1-TauD)*(1+TauI_ss+gamma*((I/K(-1))-delta));

0=min(phi,theta_ss*q*K(-1)-I);


// REST OF THE MODEL 

C=(1/h_ss)*(1-TauN_ss)*W;

p=betta*(C/C(+1))*((1-TauD(1))*Div(+1)+p(+1));

Div=(1-TauK)*(Y-W*N)-(1+TauI)*I-0.5*K(-1)*gamma*((I/K(-1))-delta)^2;

DivTau=(1-TauD)*Div; 

Div_Y=DivTau/Y;

Div_K=DivTau/K(-1);

W=(1-alp)*A*((K(-1)/N)^alp);

Rk=alp*A*(Y/K(-1));

K=(1-delta)*K(-1)+I;

Y=A*((K(-1))^alp)*(N^(1-alp));

Y=C+I;

q=betta*(C/C(1))*((1-TauD(1))*((1-TauK(1))*Rk(1)+0.5*gamma*((I(1)/K)^2-delta^2))+q(1)*(1-delta+theta_ss*phi(1)));

q_GDP=q/Y;

p_GDP=p/Y;

//TAXES AND SHOCKS

TauK=TauK_ss;

TauI=TauI_ss;

A=((A_ss)^(1-rhoA))*((A(-1))^rhoA)*exp(sd_A*epsA);

TauD=((TauD_ss)^(1-rhoTauD))*((TauD(-1))^rhoTauD)*exp(sd_TauD*epsTauD);   //AR(1) process for Div Tax Shock

theta=((theta_ss)^(1-rhotheta))*((theta(-1))^rhotheta)*exp(sd_theta*epstheta);


//Annualized Equations

YY=100*(Y/Y_ss)-100;

CC=100*(C/C_ss)-100;

II=100*(I/I_ss)-100;

KK=100*(K/K_ss)-100;

qq=100*(q/q_ss)-100;

pp=100*(p/p_ss)-100;


DivDiv=100*(Div/Div_ss)-100;

DivTau_annual=100*(DivTau/DivTau_ss)-100;

Div_Y_annual=100*(Div_Y/DivYtau)-100;

Div_K_annual=100*(Div_K/Div_K_ss)-100;

q_GDP_a=100*(q_GDP/q_GDP_ss)-100;

p_GDP_a=100*(p_GDP/p_GDP_ss)-100;

IK=I/K(-1);       //deviations in levels

IK_a=100*(IK/delta)-100;

IY=I/Y;

IY_a=100*(IY/IY_ss)-100;

CY=C/Y;

CY_a=100*(CY/CY_ss)-100;

phi_annual=100*(phi/phi_ss)-100;     //deviations in levels

N_annual=100*(N/N_ss)-100;

AA=100*(A/A_ss)-100;

theta_annual=100*(theta/theta_ss)-100;

TauD_annual=100*TauD;   //deviations in levels

end;


//// STEADY STATES DERIVED IN THE PARAMETER BLOCK

steady_state_model;

C=C_ss;
p=p_ss;
Div=Div_ss;
DivTau=DivTau_ss;
Div_Y=DivYtau;
Div_K=Div_K_ss;
W=W_ss;
Rk=Rk_ss;
K=K_ss;
N=N_ss;
I=I_ss;
q=q_ss;
Y=Y_ss;
TauK=TauK_ss;
A=A_ss;
TauD=TauD_ss;
TauI=TauI_ss;
theta=theta_ss;
phi=phi_ss;
IK=I_ss/K_ss;
IY=IY_ss;
CY=CY_ss;
q_GDP=q_GDP_ss;
p_GDP=p_GDP_ss;
YY=0;
CC=0;
II=0;
KK=0;
qq=0;
pp=0;
DivDiv=0;
Div_Y_annual=0;
Div_K_annual=0;
N_annual=0;
phi_annual=0;
TauD_annual=100*TauD_ss;
theta_annual=0;
AA=0;

end;

steady;

shocks;

var epstheta; stderr 0;     //Financial Shock not used
var epsA; stderr 0.0165;         
var epsTauD; stderr 0;                  



end;

check;

steady;


stoch_simul(order=1, periods=10000, irf=0);

